package scu.maqiang.possion;

import java.util.Arrays;

import scu.maqiang.fes.BCType;
import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3H201;
import scu.maqiang.mesh.Mesh3H20;
import scu.maqiang.mesh.Mesh3H8;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class Possion3DH201 {
	public static void main(String[] args) {
		int N = 8, label = 10;
		double ep = 1.0 / N;

		Mesh3H8 mesh0 = new Mesh3H8().cube(10, 10, 10);
		Mesh3H20 mesh3D = new Mesh3H20(mesh0);

//		mesh3D.displayNodes();
//		mesh3D.displayElements();
//		mesh3D.displayBoundarys();
//		System.exit(0);
		//
		// mesh3D.setDomainLabel(regionFunc, null, label);
		//
		FES3H201 fs = new FES3H201(mesh3D);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		//
		// GeneralCoefFunc regionCoef = (xyz, relabel, param) -> relabel == label ? 0.01b : 1.0;
		ScalarFunc ConstCoef = (xyz, relabel, param) -> 1.0e0;
		ScalarFunc[] gcf = {ConstCoef};
		//
		
		fs.assembleStiff(gcf, null, BVPType.COMMON, A);
		// //fs.assembleHeatStiff(new double[] {1.0}, BVPType.COMMON, A);
		// //System.out.println(A);
		//// System.out.println(A.makeEntrySymmetric());
		double[] RHS = new double[fs.GetNdof()];
		// //GeneralCoefFunc ConstCoef = (xyz, relabel, param) -> 1.0;
		ScalarFunc coefFunc = (xyz, relabel, param) -> 2 * xyz[0] * (1 - xyz[0]) * xyz[1] * (1 - xyz[1]) 
				 + 2 * xyz[1] * (1 - xyz[1]) * xyz[2] * (1 - xyz[2]) + 2 * xyz[2] * (1 - xyz[2]) * xyz[0] * (1 - xyz[0]);
		gcf[0] = coefFunc;
		fs.assembleSource(gcf, null, BVPType.COMMON, RHS);
		A.makeEntrySymmetric();
		fs.applyBCZero(A, RHS, Direct.X, BCType.MBN, 1, 2, 3, 4, 5, 6);
		//
		// //System.out.println(Arrays.toString(RHS));
		//
		//// System.out.println(A);
		//
		//// IterSSolver solver = new IterSSolver(1.0e-15, 1000000);
		////
		double[] x = new double[fs.GetNdof()];
		//// CSRMatrix mat = new CSRMatrix(A);
		//// solver.PCGSSOR(A, RHS, x, 1.5e0, 1);
		//// Arrays.fill(x, 0.0);
		//// solver.PCGSSOR_2(A, RHS, x, 1.5e0, 1);
		// //
		////
		//// //System.out.println(A);
		////
		//// //System.out.println(mat);
		////
		//// Arrays.fill(x, 0.0);
		//// solver.CG(mat, RHS, x, 1);
		////
		//// Arrays.fill(x, 0.0);
		//// solver.CG(A, RHS, x, 1);
		//
		Arrays.fill(x, 0.0);
		NewIterSSolver newSolver = new NewIterSSolver(A, 1.0e-10, 1000000);
		newSolver.PCGSSOR(RHS, x, 1.5e0, 2);
		//
//		Arrays.fill(x, 0.0);
//		// //newSolver.GMRES(RHS, x, 200);
//		newSolver.GMRES_SSOR(RHS, x, 1.5, 30, 1);
//		mesh3D.toTecplot("Solutions3DH20.dat", x);

        double[][] gradX = new double[3][fs.GetNdof()];
        double[][] sdX = new double[6][fs.GetNdof()];
        fs.computeFirstAndSecondDerivative(x, gradX, sdX);
        //fs.computeGradient2(x, gradX);
//        mesh2D.toTecplot("Solutions.dat", x);
//        mesh2D.toTecplot("SolutionsGradient.dat", gradX);
        System.out.println("Gradient: ");
        System.out.println("Dx min: " + MVO.min(gradX[0]) + "\tmax: " + MVO.max(gradX[0]));
        System.out.println("Dy min: " + MVO.min(gradX[1]) + "\tmax: " + MVO.max(gradX[1]));
        System.out.println("Dz min: " + MVO.min(gradX[2]) + "\tmax: " + MVO.max(gradX[2]));
        System.out.println("Second Derivative: ");
        System.out.println("Dxx min: " + MVO.min(sdX[0]) + "\tmax: " + MVO.max(sdX[0]));
        System.out.println("Dyy min: " + MVO.min(sdX[1]) + "\tmax: " + MVO.max(sdX[1]));
        System.out.println("Dzz min: " + MVO.min(sdX[2]) + "\tmax: " + MVO.max(sdX[2]));
        System.out.println("Dxy min: " + MVO.min(sdX[3]) + "\tmax: " + MVO.max(sdX[3]));
        System.out.println("Dyz min: " + MVO.min(sdX[4]) + "\tmax: " + MVO.max(sdX[4]));
        System.out.println("Dxz min: " + MVO.min(sdX[5]) + "\tmax: " + MVO.max(sdX[5]));
		// ArrayList<Integer> li = new ArrayList<Integer>();
		// li.ensureCapacity(100);
		// for(int i = 0; i < 100; i++) {
		// li.add(null);
		// }
		// li.add(1, 2);
		// System.out.println(li.toString());

		// Mesh2DT3 mesh6 = mg.square2D(2, 2, T3Type.Left);
		// mesh6.displayBoundarys();
		// mesh6.extractBoundarysFromElements();
		// mesh6.displayBoundarys();
	}
}
